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Abstract 

We investigate numerically the orbital evolution of massive extrasolar 
planets within central cavities of their parent protoplanetary discs. Assuming 
that they arrive at the inner edge of the disc due to type II migration, we 
show that they spiral further in. We find that in magnetospheric cavities more 
massive planets stop migrating at a larger distance from the edge of the disc. 
This effect may qualitatively explain the correlation between masses and 
orbital periods found for massive planets with P shorter than 5 days. 
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1 Introduction 

The standard model of giant planet formation through core-accretion and envelope 
capture assumes that a solid core is formed first by accretion of planetesimals, and 
when it becomes massive enough it begins to accrete gas from the surrounding 
nebula. It is generally believed that gas giants form exterior to the snow line, i.e. 
where the disc is cold enough for water to condense, thus increasing the amount 
of solids available for core buildup (in very massive discs the cores can reach the 
critical mass also interior to the snowline). Once a sufficiently massive planet is 
formed, it opens a gap in the disc and undergoes migration of type II; see e.g. Ida 
& Lin (2004) and references therein. Less massive planets, unable to open a gap. 
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Figure 1: The histogram of semimajor axes of Hot Jupiters. Based on data from 
Schneider (2007). 

are likely to undergo migration of type I; see e.g. Ward (1997). In either case, the 
final orbit of the planet is expected to be much tighter than the original one. 

The prime candidates for planets that have undergone extensive migration are 
objects on the tightest orbits, whose orbital radii are clustered between a ^ 0.03 
a.u. and a ^ 0.06 a.u. (Fig. [I]). A major fraction of these objects are the so-called 
Hot Jupiters (hereafter HJ), i.e. giant gaseous planets of Jupiter or Saturn type. 
Several ideas have been proposed how to stop the migration at small orbital radii, 
including presence of a central cavity in the disc, tidal interactions with the star, 
and mass loss from the planet due to the Roche lobe overflow (Trilling et al. 1998). 
The cavity may originate due to the truncation the disc by the magnetosphere of 
the star (Lin, Bodenheimer & Richardson 1996; hereafter LBR) or the magne- 
torotational instability, which accelerates the accretion flow, leading to a strong 
drop in surface density up to the radius where the temperature falls below 1500 K 
(Kuchner and Lecar 2002). 

The magneto spheric truncation hypothesis is strongly supported by the fact 
that T Tau stars have magnetic fields with intensities up to ~ 4.5 x 10^ G (Syming- 
ton et al. 2005). For typical accretion rates inferred for these objects (10^^ to 
10~ ''Mq yr~^) the disc should be truncated at a distance of a few (stellar radii) 
from the star. This prediction finds some observational support, as circumstellar 
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cavities of a few i?* are required to explain the widths of CO lines most likely 
originating in a gas on Keplerian orbits around T Tau stars (Bouvier et al. 2006 
and references therein). An obvious conclusion is that many (perhaps most) HJs 
with semimajor axes smaller than ~ 0.06 AU must have spent some time within an 
active Magnetosphere, i.e. within the central cavity of their parent protoplanetary 
disc. 

The orbital evolution of a planet within the magneto spheric cavity was consid- 
ered by LBR. In their scenario, the disc is truncated at the magneto spheric radius 
rm ~ 0.08 AU. The star corotates with the inner edge of the disc, to which it is 
magnetically locked. A planet that has migrated into the cavity experiences nega- 
tive torques from both the disc and the star. However, once it has spiralled past the 
2: 1 resonance at the orbital radius a2:i = O.GSrm, the torque from the disc becomes 
drastically weaker. The stellar torque also weakens in time due to the contraction 
of the star on its way to the main sequence (in fact, it may even reverse the sign 
if the angular momentum of the star is conserved). As a result, LBR expect the 
migration to effectively stop at ~ 0.05 AU. The same basic scenario is repeated 
in many later papers, e.g. Lin et al. (2000), Kuchner & Lecar (2002), Eisner et al. 
(2005), Santos et al. (2005) or Romanova & Lovelace (2006). However, despite 
of its popularity, it has never been verified by numerical simulations. 

LBR and their followers implicitly assume that i) the magnetosphere simply 
truncates the disc at rm, leaving it undisturbed for r > r^, and ii) the gravity of 
the matter accreting onto the star along the field lines and/or blown away from 
the star as a magnetized wind can be neglected. At a first glance, neither of these 
assumptions looks plausible; however they are quite strongly supported by sev- 
eral arguments. First, Romanova & Lovelace (2006) demonstrate that in many 
cases the magnetospheric cavity is indeed almost empty. Second, the simulations 
of Long, Romanova & Lovelace (2005; hereafter LRL) show that in many cases 
the gravity of the wind is indeed negligible (typical of their type II models is an 
outflowing corona with a density ~ 10^ times lower than the midplane density of 
the disc). Third, in some LRL models the accreting matter forms a nearly spher- 
ical thin shell with narrow polar funnels. If the shell is nearly uniform, then the 
gravitational potential it generates within the magnetosphere is nearly constant, 
and the gravitational effect of the accretion flow on the planet residing in the 
magnetosphere may be rather weak. Of course, because of numerous plasma in- 
stabilities, a smooth magnetospheric shell is an oversimplification. However, it is 
quite feasible that the time-averaged contribution to the torque on the planet from 
the fluctuating part of the density field is insignificant. 

Based on these arguments, the magnetospheric migration can be followed by 
simplified simulations in which i) the accretion disc has a sharp inner edge at 
the magnetospheric radius rm, ii) the matter that arrives at rm does not enter the 
magnetosphere but is lifted off the disc, and iii) at r < rm the net gravitational 
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effect of the magneto spheric flow is neglected. Following the implicit assumption 
made by LBR, we also assume that beyond the disc is practically undisturbed 
by magnetic fields, i.e. that non-Keplerian effects in its dynamics are dominated 
by the gravity of the planet rather than Lorentz force. The latter assumption, whose 
validity we check a posteriori in Sect. 12. 2[ allows us to calculate the torque from 
the disc on the planet with the help of standard two-dimensional simulations. 

The simulations are presented in Sect. 2. In Sect. 3 we discuss their results, and 
speculate about the origin of the semimajor axes distribution shown in Fig.[IJ The 
paper is supplemented by an appendix containing a report on the code calibration 
test. 



2 Numerical simulations 

We employ the FLASH code, which operates on a block-structured grid (Fryxell 
et al. 2000). Within each block the equations of hydrodynamics are solved using 
a conservative third-order Eulerian scheme (Colella & Woodward 1984). FLASH 
can build a hierarchy of progressively finer grids by "halving" the blocks in each 
coordinate whenever refinement is required. For the present simulations, however, 
a uniform (nonrefined) grid proved to be sufficient. 



2.1 Problem setup 

We consider a system composed of a central star with mass (in the simulations 
= IMq), a planet with mass Mp, and a thin non- self gravitating disc truncated 
at a distance from the star. The disc is locally isothermal, and the vertically 
integrated pressure p is related to the surface density S through 

P = ^cl, (1) 

with the local sound speed given by 

Cs = H^K , (2) 

where ilx = a/ GM/r^ is the Keplerian angular velocity, and H is the half- 
thickness of the disc. The aspect ratio h = H/r is assumed constant over r and 
equal to 0.05. 

We do not introduce any explicit magnetic field. Following the arguments pre- 
sented in Sect. [H we simply assume that for r < the accretion flow does not 
influence the planet in any way. With such an assumption, in a coordinate system 
centered on the star the orbital evolution of the planet is described by the standard 
equation 

d^_G(M^_ (3, 

at a-^ 
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where a is the position vector of the planet, and $d is the gravitational potential 
of the disc. The formula for the potential, 



S(r) 



■dr + G 



r — a 



— —a ■ rar , 



(4) 



contains the indirect term accounting for the non-inertiality of the coordinate sys- 
tem, and the integration in dH) is performed over the surface of the disc. 
The magneto spheric radius is given by the standard formula 



GM^M2 



1/7 



(5) 



where 5*, R-^, M^, M and r] are, respectively, field strength at the surface of the 
star, stellar radius, stellar mass, accretion rate, and a dimensionless factor of order 
unity. For a star with an aligned dipole field 0.5 < 77 < 1.0; see e.g. Lai (1999) and 
references therein. Normalized to the standard parameters of T Tau stars, equation 
^ reads 
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We take = 12Rq = 0.056 AU, which agrees with the distance range suggested 
by the distribution shown in Fig. [U (note also that for the standard values of B^,, 
and M the adopted value of rm corresponds to an entirely reasonable range of 
stellar radii 2Rq < R^ < 3Rq). 



2.2 Results 

The simulations were performed on a polar grid of Ur = 512 and n,^ = 128 
points, extending from rin = to rout = 0.356 AU. We checked that doubling 
Ur and changes the value of the disc torque on the planet by less than 2%, 
which proved that our basic resolution was sufficient for the present problem. At 
Tin a free outflow boundary condition was applied, while at the outer edge of the 
disc we implemented the prescription proposed by de Val-Borro et al. (2006) to 
damp reflections from the boundary of the grid. Specifically, after each time-step 
the equation 

was solved, where X stands for S or velocity component, Xq is the initial profile 
of a particular variable, and R{r) is a parabolic ramp function which decreases 
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Figure 2: The migration of planets in a disc with Md = lOMj. Solid: Mp = 0.5Mj. 
Dashed: Mp = IMj. 

smoothly from 1 at r = rout to at r = 0.84rout- The default FLASH artificial 
viscosity with Cvisc = 0.1 was used to stabilize the solution. 

Initially the disc is uniform, and a = 0.062 AU. To relax the model, we keep 
the planet on the initial orbit until it carves a quasistationary gap in the disc, ex- 
terior to which a spiral density pattern emerges. When the relaxation procedure 
is completed (usually after ~ 100 orbital periods of the planet), the time-counter 
is reset to 0, and the planet begins to evolve according to equation (|3]) which we 
integrate with a fourth-order Runge-Kutta scheme. We performed simulations for 
two values of the mass of the planet (Mp = 0.5Mj and Mp = IMj, where Mj is 
the mass of Jupiter), and two values of the mass of the disc contained in the grid 
(Md = lOMj and Md = lOOMj). The latter were chosen unrealistically large in 
order to speed the evolution up, so that noticeable effects could be produced in 
a reasonable CPU time. Within the adopted scenario the gravitational torque T 
from the disc on the planet should scale linearly with Md, and this is what we 
observe in our simulations. Thus, the fact that our discs are far too massive does 
not influence the conclusions of the paper. 

As soon as the planet has been released from its initial orbit, it begins to mi- 
grate into the magneto sphere. The migration process is illustrated in Figs. [2] and [3] 
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Figure 3: As in Fig. O but in a disc with Md = lOOMj 

for discs with Md = lOMj and Md = lOOMj, respectively. As expected, a planet 
with given Mp migrates faster in more massive discs. Figs. [21 and [3] also demon- 
strate that in a disc with given Md the migration rate is faster for more massive 
planets. 

Finally, in Fig.|4]we plot |T| as a function of r for both values of Mp (note that 
after an initial adjustment the torque scales with Mp). In the following we argue 
that for realistic values of A/d,p the gravitational torque is comparable to or weaker 
than the magnetic torque Tm on the disc from the star. In the following, T and Tm 
will stand for absolute values of gravitational and magnetic torque, respectively. 

The centrally peaked minimum mass solar nebula (MMSN) with the surface 
density profile 

Au) 

(Ruden 1999), and the more gently peaking disc simulated by Nelson et al. (2000), 
with the profile 

E2 = 10=^gcm-2(^)"'^' (9) 

yield, correspondingly, Md ~ 1.5Mj and Md ~ 0.03Mj. For further estimates 
we take the geometrical mean of these two values, i.e. Md = 0.2Mj. The standard 
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Figure 4: The absolute value of the disc torque on the planet as a function of the 
distance from the central star. Triangles: Mp = lAfj. Circles: Afp = 0.5Afj. Thin 
line: torques for Mp = 0.5Mj multiplied by 4. Vertical axis: values obtained from 
the simulation with = lOOMj, scaled to Ma = 0.2 Mj. 



formula for 



m 



= ^ (10) 



(Lai 1999) contains weakly constrained parameters B^, and i?^. However, using 
dS]) with = 1 we see that the characteristic time scale 

_ Jp(r^) _ Mp^GM^r^rj _ Mp 

" B^^Rt "m ^''^ 

(where Jp(rni) is the orbital angular momentum of the planet at the edge of the 
magneto sphere), does not depend on 5^ or i?^. For Mp = IMj and M = 1O~'^M0 
yr^^ equation (fTTI) yields Tm = 10^ yr. This is the time scale on which a planet 
placed at would migrate if the torque from the disc was equal to T^. We shall 
compare it with the timescale of migration due to the gravitational torque, Tg = 

The orbital angular momentum of a IMj planet placed at = 0.056 AU is 
equal to 1.5 x lO^^MgAU^yr^^ From Fig.Hwe have T{r,n) = 7.5x lO^^MgAU^yr" 
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yielding Tg = 2 x 10^ yr. Thus, T = 5Tm, and the assumptions of Sect. [T] are sat- 
isfied. However, already when the planet has migrated to a = 0.9r^ = 0.05 AU, 
T drops to lO^'^MoAU^yr^^, and Tg increases to 1.5 x 10^ yr. This means that 
is now larger than T, i.e. the inner edge of the disc is more strongly influenced by 
MHD effects rather than gravity. For the 0.5Afj planet the situation is even worse, 
since T ^ already at r^. As the disc dominated by chaotic MHD effects can- 
not exert a consistent negative torque on the planet, we conclude that planets with 
Mp < 0.5Afj are unable to spiral into the magnetosphere, whereas Jupiter-size 
planets stop spiralling at a distance ~ 0.006 AU from the inner edge of the disc. 

3 Discussion 

Based on the data from the beginning of September 2004, Mazeh, Zucker and Pont 
(2005); hereafter MZP, and Gaudi et al. (2005) found that among 6 then-known 
transiting HJs more massive planets tended to have shorter periods. The same 
effect was visible in the whole population of 17 then-known planets with orbital 
periods shorter than 5 days (hereafter: SPP), albeit with large scatter. Presently 
(end of May 2007) 46 SPP are known (47*'' object with P < 5^^ is a brown dwarf 
HD41004 B b), 20 of which are transiting (Schneider 2007), and it is interesting 
to see whether the original correlation still persists. 

Following MZP, we exclude planets in known binary systems (r Boo b and 
HD 188753A b) as weU as "Hot Neptunes" (GJ 436 b, 55 Cnc e , Gliese 876 d, GJ 
674 b and HD 219828 b), which, as MZP write, "probably are of a different nature 
and have a different formation and evolutionary history". We also exclude three 
transiting objects: SWEEPS-04, whose RV amplitude is insignificant, SWEEPS- 
11, which seems to be a brown dwarf rather than a planet (Sahu et al. 2006), and 
XO-3, whose mass is at the lower limit of the brown-dwarf range, and whose or- 
bit is highly eccentric. Amazingly, the slope of the linear fit to the Mp(P) relation 
found for the 14 transiting planets with reliably determined masses is nearly the 
same as the one found by MZP (Fig 15]). Moreover, the correlation is also visible 
for the whole sample of SPP, and the slope also nearly the same. MZP could not 
provide any explanation for the observed effect, and suggested that such an expla- 
nation should be worked out when the Mp — P correlation is better established. As 
the new data confirm their findings, the effect seems to be ripe for consideration. 

Because masses of stars orbited by SPP are rather tightly clustered around 1 
Mq, whereas masses of SPP are at least several hundred times smaller, periods 
of SPP are primarily determined by their orbital radii. As a result, the Mp — P 
correlation implies a similar connection between masses and orbital radii, which 
is indeed observed (Fig. [6l). We have shown that the "classical" scenario of the 
orbital evolution inside the magnetosphere, in which planets stop migrating at 
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Figure 5: The mass-period relation for planets with periods shorter than 5 days. 
Based on data from Schneider (2007). Squares: transiting planets. Triangles: plan- 
ets for which spectroscopic data are only available (as in MZP, their lower mass 
limits are divided by 7r/4, the expected value of sini). Dotted, dashed and solid 
line: fits to, correspondingly, 6 transiting planets of MZP, 17 transiting planets 
with reliably determined masses known in January 2007, and all SPP known at 
the end of May 2007. 



~ 0.63rjn, is not true. Moreover, even if it were correct, it would not be able 
to account for the observed effects. To see it, imagine an evolving population of 
planets represented by points scattered on the (a, Mp) plane. Initially, every planet 
resides at its r^. We do not know the distribution of rm, but it seems reasonable 
to assume that it is flat between some rjn.min < 0.056 AU and ^ax ~ 0.056 
AU. If this is so, then 1) at t = there's no correlation between Mp and a, and 
2) such correlation cannot emerge for any t > 0, because every point is displaced 
horizontally by an amount which depends on r^, but does not depend on Mp. 

If the cavities are not of magnetospheric origin, but originate from increased 
viscosity and accelerated accretion flow (Kuchner and Lecar 2002), we encounter 
another problem. Although in this case the disc is not perturbed by the magnetic 
field of the star, its inner part may be strongly unstable, causing the edge of the 
cavity to irregularly or semi-regularly change its location (Wiinsch et al. 2006). 
As a result, any correletion between Mp and a or P is highly unlikely. 
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The results of our simulations imply however that more massive planets tend 
to stop deeper in the magnetosphere, and in the following we argue that this ef- 
fect can at least qualitatively account for the observed correlations. If the disc 
conforms to the simple a-model, then the velocity of the accretion flow is given 
by 

Vj- = a/i^ri^K, (12) 
and, provided that a and h do not change, the accretion rate 

M = 27ra/iV2fiKS (13) 

scales proportionally to S. This means that by varying M at a constant r^a we do 
not change the timescale ratio t^/t^, i.e. the conclusions of Sect. [21 remain valid 
when the surface density is decreased or increased at the inner edge of the disc. 

As the planet recedes from the inner edge of the disc, its tidal effect on the 
disc approaches the linear regime in which the density perturbation AS is small 
compared to the unperturbed density Sq. Eventually, T should converge to the 
linear torque which scales proportionally to aA/p (Ward 1997), and Fig. |4] shows 
that it is indeed the case. Based on a particle approximation, the same scaling was 
found by Lin and Papaloizou (1978) for the case of a small-mass component of 
a binary with an external (i.e. circumbinary) disc. (Strictly speaking, the formula 
given by Ward concerns the net torque on the planet, but the torque components 
originating exterior and interior to the orbit scale in the same way). After some 
algebra we see that timescale ratios of planets with masses M^i and Afp which 
begin to spiral into cavities whose edges are located, correspondingly, at and 
rjn, are related by the formula 



Tgi _ ( Mp\ pF^ Tg 



(14) 



pi / V ™l ™ 



Since we have found that ~ O.SAfj-planets are hardly able to detach from 
the edge of the cavity, we may assume based on Fig. [6] that the radii of the 
smallest cavities were not much different from ~ 0.036 AU. Let us consider 
two identical planets, the first one spiralling into our "standard" magnetosphere 
with rm = 0.056 AU, and the second one into the "compact" magnetosphere with 
^'mi = 0.036 AU. According to equation (fT4)) . for the second planet the ratio Tg/ 
is just by ~ 25% larger than for the first one. As a result, both planets reach a sim- 
ilar stopping distance from the edge of the disc, and to a first approximation we 
may assume that the function (is (Afp, does not depend on r^. 

According to Marcy et al. (2005), the distribution of planetary masses is af- 
fected very little by the unknown inclination of orbits, and based on available data, 
one may adopt 

dN{Mp) M-^dM^. (15) 
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Figure 6: The mass - semimajor axis relation for the same planets as in Fig. |51 
Solid, dotted, and dashed line: fits to, correspondingly, all SPP known at the end 
of May 2007, final (a, Mp) distribution for dN{Mp) ~ M-^dMp, final (a, Mp) 
distribution for dN{Mp) ~ M'^dMp. 



As we have mentioned before, the distribution of magneto spheric radii is un- 
known, but since there is no obvious reason to favour any value of rm we may 
assume that all values between r\ai and equally probable. With these as- 

sumptions, the joint distribution of Mp and rm is given by 

F(Mp,r^) =ArMp-ixi(Mp)x2M, (16) 

where Xii^) = 1 for 0.5 < x < 2.5 and otherwise (2.5 Mj is the estimated 
mass of the most massive SPP in Fig.|6l), X2ix) = 1 for 0.036 < x < 0.056 and 
otherwise, whereas 

(POO POO \ —1 

/ / m^^ xii'^TT') X2{x) dx dm ] . (17) 

The ratio t^/t^ decreases rather sharply with the increasing mass of the planet, 
but in order to remain on the conservative side we may adopt that the stopping 
distance scales linearly with Mp, according to the formula 
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which in agreement with our results yields dg = AU for Mp = 0.5Mj and 
4 = 0.006 AU for Mp = IMj. 

Let us again imagine an evolving population of planets represented by points 
scattered on the (a, Mp) plane. Every planet begins to spiral down from its r^, and 
stops at cis = — ds, with ds depending primarily on the planet's mass. Thus, the 
final distribution of points is skewed toward smaller a, i.e. a correlation betwen 
Mp and a is generated. The least-square linear fit to the final distribution 

Mp = a * a + /3 (19) 

is given by the standard formulae 

/o°° /o°° ("^ ^ ^p) * [(^ ~ ds{m)) — (fm — ds)] F{m, x) dx dm 
/o°° JT [(^ ~ ds{m)) - (fm - ds)Y F{m, x) dx dm 
13 = Mp-a* {f^-ds) , (20) 

and is shown in Fig.[6l For comparison, analogous fit obtained for a steeper mass 
distribution dN(Mp) ~ M~'^dMp is also plotted. Given the simplicity of our 
approach, the qualitative agreement of theoretical fits with the observational one 
looks quite encouraging. 

So far we have neglected tidal interactions between the planet and the star. 
Following LBR and recent theoretical as well as observational evidence (Bouvier 
et al. 2006; Herbst & Mundt 2005, and references therein) we assume that T Tauri 
stars are rotationally locked to inner edges of their accretion discs. Since a locked 
star rotates at fi^ ^ ^(^m) < the resulting stellar torque onto the planet is 
negative, causing orbital radius of the planet to evolve according to the formula 

a{tr = aiOr - Y^^^^ ' (21) 

see e.g. Piitzold et al. (2004). For M^ = Mq we have 

x{tr = xr- 3.2x10-'^^ i^^^t, (22) 

where x = a/ Rq, t is expressed in years, and a value of 1.2 x 10^ was adopted 
for Qi,/k2 (Patzold et al. 2004). Under the influence of stellar tides alone a 1 Mj 
planet departing from the edge of the magnetosphere at = 0.056 AU would 
spiral down to the surface of a 3i?0 star within ~ 1.3 x 10^° yr, whereas for the 
same planet departing from r^i = 0.036 AU the spiralling would take ~ 7.5 x 10^ 
yr. In fact, the spiralling time would be even longer because the star contracts, 
causing the stellar torque on the planet to decrease. According to AUain (1998), 
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within ~ 10^ yr R^, can easily shrink by a factor of 5 if the star does not accrete, 
and when accretion is allowed for the contraction proceeds even more quickly 
(Siess et al. 1999). As a result, it is highly unlikely for stellar tides to significantly 
influence the orbital evolution of the planet within the disc lifetime. 

On the other hand, equation (|2T]) predicts a faster evolution for more massive 
planets, and assuming that Vl^ remains smaller than one might expect that Fig. 
[5]is just a snapshot of a distribution undergoing slow secular changes due to stellar 
tides. However, at the end of the pre-Main Sequence phase the radius of the star is 
equal to ~ Rq, and the stellar torque on the planet strongly decreases. As a result, 
the spiraling timescale of a 1 Mj planet from r^i = 0.036 AU lengthens to well 
above 10^° yr. 

In the final note we would like to stress that we do not claim that the prob- 
lem of correlation between SPP masses and periods has been solved. We merely 
report a mechanism which may qualitatively account for the observed trend. An 
alternative mechanism was disussed by Faber, Rasio & Willems (2005) and Rasio 
& Ford (2006), who argue that SPP are likely to have originated from the tidal 
capture of planets on originally highly eccentric orbits. It is also possible that both 
mechanisms have been at work; provided, of course, that the correlation is real 
and will be confirmed by future observations. 
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Appendix: Code testing 

Before the version of FLASH adapted to the problem of disc - planet interac- 
tion was used for the proper simulations, we had tested it for compatibility with 
analytical expressions for the disc torque on the planet obtained by Ward (1997). 

Consider a thin disc with a surface density a{r, 0), orbiting a central star of 
mass M^. The mass of the disc is minute compared to M^, and its self-gravity is 
negligible, causing orbital velocity and epicyclic frequency of the disc matter, n 
and K, to be nearly Keplerian. Assume that the temperature distribution in the disc 
does not depend on time, with the local sound speed Cg fixed according to The 
disc is perturbed by a planet of mass Mp = /xM^, which moves around the star on 



16 S. Starczewski, A. J. Gawryszczak, R. Wiinsch & M. Rozyczka 



a circular orbit of radius a. As the perturbations are nonaxisymmetric, the planet 
is subject to a gravitational torque T from the disc. According to Ward (1997), the 
torque density dT/dr is given by 



dT _ 2ij\aa^){anpfmf^^ fn 
d7 ~ ^ r(l + 4^2) 



(23) 



where Vlp is the angular orbital velocity of the planet, ^ = m^Cs/rK, and e is equal 
to +1 (—1) for torques originating in the disc interior (exterior) to the planet's 
orbit. The function mr is given by the formula 



rrir 



(n - i]p)2 



c2/r2 



1 



-1/2 



(24) 



where Q,'^ 
formula 



fi2 



cl/r'^ and ^ were used; whereas the function - by the 
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1 



^1/2 
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+ 2Vl + e hr,2{x 



(25) 



where &™2(^) Laplace coefficient with argument x = r/a. 



To calculate the torque density from the numerical models, we integrate con- 
tributions from all grid cells forming a thin annulus at a distance r from the star : 



dT 
dr 



where 



G 



d 



2tt 



Mp r sin( 



r cos 



(26) 



(27) 



is the distance between the planet and the point (r, </>) on the annulus; (pp stands for 
the positional angle of the planet, tr is the Roche radius, and e is a dimensionless 
softening parameter. 

For the test runs we adopted Mp = 10 Me (Earth mass) and a = r^. The 
calculations were performed on a polar grid (r, 0) extending radially from 0.03 to 
0.18 AU. The disc had a constant surface density a = 2.28 x 10^ gcm^^, so that 
the corresponding mass of the disc matter contained within the grid was equal 
to 2.5 Mj (Jupiter mass). The softening parameter e was set to 0.6. No explicit 
viscosity was used, and no special treatment was applied to the matter within 
the Roche lobe of the planet. The disc gas could flow freely through the radial 
boundaries of the grid. 

We ran two test simulations with different resolutions. The results are shown 
in Fig.|71 from which it is evident that on the finer grid (512 x 256 points) the code 
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Figure 7: Torque densities obtained numerically for Ur x n^j, = 256 x 128 points 
(dotted) and 512x256 points (dash-dotted), compared with the analytical so- 
lution (solid). Vertical lines: location of the planet (solid) and the torque cut- 
toff limits at the Roche radius (dashed). Tq is the normalized torque defined as 
To = 7t/i^(7a^(ar2p)^(a//i)^. The plotted values are averages from 100 time steps 
(roughly one radian of the orbit). 



is able to calculate reasonably accurate torques from the parts of the disc that are 
more distant than ~ 0.005 AU or ~ 10 grid cells from the planet. Since that area 
contributes less than 20% of the total torque from the disc, it may seem rather 
discouraging. The obvious conclusion is that if the planet is not able to open a 
gap, a really very high resolution is necessary to follow its migration through the 
disc. Obviously, for gap-opening planets the resolution requirements are much 
less stringent. 



